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This paper treats the basics of thermoacoustic engines. The set of differential equations, 
which describes the dynamics of the individual components, is condensed in a single 
high-order differential equation which determines the time dependence of all dynamic 
variables. From this relation analytical expressions are obtained for the damping 
coefficient, the oscillation frequency, and the onset temperature that allows stable 
oscillations. Also transient effects are discussed based on numerical integration of the 
dynamic equations. 

© 2009 Elsevier Ltd. All rights reserved. 


1. Introduction 

This paper treats the basics of thermoacoustic engines. The formalism that is presented here is general and can be 
applied to a large variety of thermal engines. The basic idea is that the set of differential equations, which describes the 
dynamics of the individual components, is condensed in a single high-order differential equation which is the heart of the 
dynamics of the entire system. It determines the time dependence of all dynamic variables. The so-called travel ling-wave 
thermoacoustic Stirling engine, described in the paper by Backhaus and Swift [1 ], will be taken as the model system. This is 
a very important type of thermoacoustic engine which is studied in many laboratories around the world. Fig. 1 is a 
schematic drawing of this engine. It consists of a long resonator tube and a loop which contains a regenerator, several heat 
exchangers, a compliance (c), a connecting tube (d), a pulse tube (t), and a section, with a smaller diameter, called the 
inertance. One of the characteristic features of this engine is that the gas starts to oscillate back and forth in the machine 
“spontaneously” if the middle heat exchanger, with temperature T t , is heated to a temperature which is sufficiently high 
above ambient temperature T a . This effect is the main topic of this paper. 

After explaining the model, the various differential equations will be derived, taking the temperature T t as a 
fixed input parameter. This results in a single differential equation for the pressure p t in the pulse tube. It will turn 
out to be a fourth-order differential equation with real constants. It will subsequently be solved analytically and its 
properties will be discussed. Expressions for the damping coefficient and conditions for the onset and stability of 
oscillations are derived as well as the relations for the oscillation frequency. In the second part of the paper the heating 
power Qt is fixed and T t is treated as a slowly varying function of time. For this case the dynamic equations 
are solved numerically. It reveals that the state of steady oscillations is reached after an interesting transient period 
which may show bursts of high-amplitude oscillations. So far this effect has been poorly investigated experimentally. In the 
state of steady oscillations the results of the analytical treatment and the numerical calculations are in excellent 
agreement. 
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Nomenclature 

Zr 

geometrical flow impedance of the regenerator 




(Eq. (27)) (nr 3 ) 

a 

notation (Eqs. (9), (38)-(41)) 



a, b, cj, k, l, m, n notations (Eqs. (58) and (59)) 

Greek symbols 

A 

area (m 2 ) 



c 

velocity of sound (m/s) 

oc 

dimensionless damping parameter 

q 

dimensionless parameters (Eq. (75)) 

y 

specific heat ratio 

C 

flow conductance (m 3 /Pa) 

n 

viscosity (sPa) 

Ch 

heat capacity of heat exchanger (J/K) 

Ka 

effective thermal conductivity (W/Km) 

D 

diameter (m) 

V 

frequency (Hz) 

f 

volume rate of change (m 3 /s) 

T 

dimensionless temperature 

* 


CD 

angular frequency (rad/s) 

H 

enthalpy flow rate (W) 

CD X 

dimensionless angular frequency 

L 

length (m) 



M 

mass (kg) 



P 

pressure (Pa) 

Subscripts 

Q 

heating power (W) 



s 

rate parameter (Eq. (42)) (s -2 ) 

0 

time average, characteristic 

t 

time (s) 

a 

ambient 

T 

temperature (K) 

ac 

acoustic 

V 

velocity (m/s) 

c 

critical 

V 

volume (m 3 ) 

c,d,f,R,t spaces in the system 

* 

V 

volume flow rate (m 3 /s) 

e 

effective 

w 

compliance factor (m 3 /Pa) 

0 

orifice 

X 

position (m), dimensionless time 

P 

(time) period 

Zr 

specific flow impedance (m -2 ) 

r 

regenerator, resonance 

Z 

mathematical parameter 

ref 

reference value 


2. Modeling 

The dimensions of the components in the loop are all sufficiently smaller than the wavelength of the oscillations so it 
can be modeled as in Fig. 2. The loop is decomposed in compartments c, t, and d (in which the pressures are homogeneous), 
an inertance, and a regenerator. The inertance is supposed to behave as a solid piston with a mass M, equal to the mass of 
the gas contained in the inertance part. The orifice, together with the buffer volume b , represent a load. We assume that 
there is dissipation only in the orifice and in the regenerator. 

All volume flows in the loop can be determined from the local pressure and its time derivative except the volume flow 
V R . The resonator is so long that there are fundamental delays in it due to the finite velocity of sound. It is a challenge for 
future analysis to model the resonator tube in acoustic terms. In this paper the resonator will be modeled as a cylinder with 
volume V R (average length L R0 and cross sectional area A R ) of a space R , closed at the right by a piston with mass M R . By 
modeling the system in this way we really discuss a kind of free-piston Stirling engine. It is assumed that the pressure 


inertance 



Fig. 1 . Schematic drawing the travelling-wave thermoacoustic Stirling engine. The symbols are explained in the text. 
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Fig. 2. Schematic drawing of the thermoacoustic engine. Compared to Fig. 1 the resonator tube is replaced by a volume V R , closed by a piston with mass 
M r . An orifice with flow conductance C 0 , backed by a buffer volume b (which is assumed to be very large), is included to take into account dissipation in 
the system outside the regenerator. 


inside V R is only a function of time. The value of M R is taken equal to the mass of the gas in V R , so M R = PqL ro A r with p 0 
the average density of the gas. The lower index 0 is used to indicate time-averaged values. The resonance frequency co R of 
the mass-gas-spring system formed by M R and V R is 



IP o 
PqLro 



The basic angular resonance frequency co ac for an “open” tube (quarter wavelength), as in the case of Ref. [1], is given by 


tOac^ac 1 

-= n 71 

c 2 

with c the velocity of sound. If we take the average length L R0 of the cylinder in such a way that co R 
Poc 2 = vpo’ 



co ac we have, using 




3. Dynamic equations 


3.1. Governing equations 


The properties of the resonator get index R and of the inertance index i. The pressures in the spaces t, d, and R are all 
equal to the pressure p t in the pulse tube t. The pressure in c is p c . The equation of motion of the mass of the inertance is 



d 2 Xj 

dt 2 


A (Pf - Pc) 



with Xj the position of the mass M i} counted positive in the upward direction. The volume of the space below M, (space d ) is 
given by 


^d ~ ^dO + A*i- 

With Eq. (5) relation (4) can be written in terms of the volume V d 

d 2 V d Af 

dt 2 Mj Pr ’ 

where we wrote p r for the pressure drop over the regenerator (and over the inertance) 


( 5 ) 

( 6 ) 


Pr = Pt ~ Pc- 


( 7 ) 




























































A.T.A.M. de Waele / Journal of Sound and Vibration 325 (2009) 974-988 


977 


We will also use the notations 


Spt =Pt-Po 


( 8 ) 


and 


4 , 4 

a R = —p~ and a,- = T f- 
m r M i 

for convenience. 

The acceleration of the mass M R of the resonator is given by 

Mr - A R 3p t . 


(9) 


Since 


we have 


V R = V R0 +A R X R 


( 10 ) 


( 11 ) 


d z V 


R 


dV 


= (d R Sp t . 


( 12 ) 


The volume flow through the orifice is given by 


^b — Co^Pt* 


(13) 


with C 0 the flow conductance of the orifice. 

In Fig. 3 three situations are depicted: one where gas flows into and out of a control volume with volume V and pressure 
p(t), one where gas flows into a control volume with a moving piston, and one where gas flows out of a container through a 
valve with flow conductance C. For these three situations we have three analogous relations. For the case of Fig. 3a 


* * V dp 

Vi = v 2 + p 


yp d t 


(14) 


The nice thing about Eq. (14) is that it holds for an adiabatic ideal gas in the container, even if the temperature is not 
homogeneous. It holds, in particular, for the pulse tube where the temperature on one side is T t and on the other side T a . 
For the case of Fig. 3b holds 


= vA + 


V dp 
yp dt 


(15) 


and for the case of Fig. 3c 


0 = C(p - p 0 ) + 


V dp 


(16) 


yp dt' 

In the analysis it will be assumed that the pressure variations are much smaller than the average pressure. In that case V/yp 
can be replaced by their average values V 0 /yp 0 . Each of the spaces R , c, d, and t in the system of Fig. 2 can be characterized 
by the parameters 

ypo 


Wi = 

1 V,- 


with i = P, c, d, t. 


(17) 


i0 


a 



pioy 




p{t\V 


* 

V 


c 


Fig. 3. Gas flows with varying pressure: (a) inlet and outlet; (b) inlet and piston; and (c) fixed volume with valve with flow conductance C. 
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With Eqs. (14)-(16) we can express the various volume flows, as defined in Fig. 2. With Eq. (14), applied to the pulse 
tube t, 



* * 1 d n t 

v h = v t + ,7 . 

(18) 

With Eq. (15) and (5), applied to volume d, 

* _AV d 1 d p t 
d dt w d dt 

(19) 

and, using Eq. (11) for R, 

,*/_ _ dV R , 1 d Pt 

R d t w R dt 

(20) 

and also for c 

* _dV d 1 d p c 
c dt w c dt ' 

(21) 

Mass conservation gives 

* * * * 

v t = v b + v d + V R . 

(22) 

For simplicity we assume that the void volume in the regenerator is 
conservation over the regenerator gives 

zero (zero compliance). In that case mass 


* * 

V h = T t V c 

(23) 

with the reduced hot-end temperature 

T - Tt 

Tf T' ’ 

1 a 

(24) 

In a linear approximation the volume flow, entering the regenerator, is proportional to the pressure drop p R and can be 
written as 


V c — —CrPr • 

(25) 

The flow conductance C r is written as 

c r = 1 

da^r 

(26) 

with r/ a the viscosity at room temperature and 

^ Jr 

N ^ 

II 

N 

(27) 


where z r is the specific flow resistance of the regenerator and L r and A r are the length and cross sectional area of the 
regenerator, respectively. In general C r depends on the temperature distribution, etc. in the regenerator. More elaborate 
treatments of the regenerator can be found in Refs. [1,2]. For more detailed treatments numerical models can be used such 
as Sage [3], Regen3.2 [4], or DeltaE [5]. 


3.2. One dynamic equation 


In the previous section the dynamic equations are given for each of the components of the engine. In this section these 
relations will be combined to the extreme, resulting in a single differential equation for one parameter (in fact 3p t ) which 
determines the over-all dynamic behavior. The procedure is to eliminate subsequently all but one variable. 

The volume flow V h in Eq. (18) can be eliminated with Eq. (23) 


* 

and, with Eq. (25), V c can be eliminated 



1 dSp t 
w f dt 


(28) 


T tCriPc 


Pt ) = V t + 


1 ASp t 
w t dt 


(29) 














A.T.A.M. de Waele / Journal of Sound and Vibration 325 (2009) 974-988 


979 


Eqs. (13), (19), (20), and (29) in (22) give 


TtCrPl+ „, 5 |. + «p, + !!& + -L^ + 5 & + -L^,o. 


With the notation 


W e = 


d t w d d t d t w R d t 


ypo 


Vt + ^dO + 


we get from Eq. (30) 


ddPt 
d t 


dV d dV r 


Eqs. (21), (25), and (7) give 


— TfWgCrP r + WeCodpt + Wg ——K We 


r n - dv d , 1 dS Pt 1 d f*R 
rPr dt w c dt w c dt ' 


dt 


With Eq. (32) 


dp, 


dV 


— (Wc + We) “T7~ + (TfWeCV + WcCj-)p r + WeCoSpt + We 


dV 


K 


dt v L 1 dt 1 v 1 " ' 1 L 1 e u 1 e dt ' 

We eliminate V R and V d by differentiating Eqs. (34) and (32) once. Substitution of Eq. (6) and (12) gives 

n 

dSu d v du 

-w e C 0 -fff - w e a R Sp t = —+ (T(We + W c )C r -gp + (W c + W e )ajP r 


and 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


d 2 Sp t _ d 3p t s „ dp r 

+ w e C 0 + w e a R Sp t -t t w e C r - w e a,p r . 


dt" 


dt 


(36) 


Eqs. (35) and (36) form a set of two linear second-order differential equations with Sp t and p r as variables. In Appendix A it 
is shown that elimination of p r leads to a single homogeneous fourth-order differential equation in Sp t 


d 4 Sp t d 3 Sp t d 2 Sp t d 3p t s _ 
+ a 3 -1" a 2 ~~77o -^ a l + a 0^Pt = 0 


dt 4 


dt 3 


dt 2 


dt 


with 


fl 3 — WeCo + WeCp + TfWeCp, 

(22 — 0 R We H - a,Wc + fl,We + WcWeCpCo, 

= WcWe(CrQ R + Cofl,), 
a 0 = w c Wea R a,-. 

Note that all factors are positive and that only a 3 depends on the dimensionless temperature T f . 


(37) 


(38) 

(39) 

(40) 

(41) 


3.3. Stable oscillations 


Eq. (37) generally results in solutions which are oscillating in time. This is shown in more detail in Appendix A. The 
function s, defined in Eq. (84) as 


s 


CL\ 

a 3 


+ a o 


03 

cq 




(42) 


determines the growth or decay of oscillations. As shown in Appendix A the function s is proportional to the damping 
coefficient a with a positive proportionality constant. If s>0 the oscillations grow, if s<0 the oscillations are damped. In 
any case s< 0 if T t = 1, so the oscillations are always damped in an isothermal system (as it should). 

For steady oscillations s = 0 so 


0 = a 0 



«2 



(43) 


or 


a 3 a 2 ± \J a 2 - 4a o 

2a 0 


( 44 ) 
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and also (Eq. (86)) 



a-[ 

os' 


So the oscillation frequency is given by 


_L /£T 

2n\ja 3 ' 


( 45 ) 


(46) 


In Eq. (44) only the solution with the minus sign gives a positive dimensionless temperature T t . Substitution of Eqs. 
(38)-(41) in Eq. (44) gives the relation for the temperature T c at which the oscillations are stable 



1 

2 W e 



Cq We 
Cj- We 


(47) 


If the system is originally at room temperature and a sufficiently large heating power is applied, the oscillations start if the 
temperature has reached the critical value of T c . Therefore this temperature will be called the onset temperature. 
Substituting t c in Eq. (38) and (45) gives the angular frequency of the oscillations. 

Substituting Eqs. (39) and (41) in Eqs. (47) and (38) and (40) in Eq. (46) leads to rather complicated expressions which 
show how the various system parameters affect the onset temperature and the oscillation frequency. If C 0 = 0 
(no dissipation outside the regenerator) t c is given by 


ay \/ai - 4a n - 2w c a,- 

Tc =—- ifc °= a ^ 

It is independent of the regenerator conductance C r . This somewhat surprising result is due to the fact that, on one hand, 
the flow in the regenerator gives dissipation, but, on the other hand, the pressure difference over the regenerator is the 
driving force for the oscillations. In the limit of small w e (or large V R0 ) and with C 0 = 0 a second-order series expansion of 
the numerator in w e shows that the onset temperature is given by 


Tc _ w e a R 
T a ~ W c Oj 


(49) 


and the corresponding frequency by 


oj 2 = w e a R = oj 2 r . (50) 

So, in this limiting case, the system oscillates at the resonance frequency of the resonator tube. With the values of Table 3 
Eq. (49) would lead to values of t c close to one while in practice t c is more on the order of two. This shows that, in reality, 
dissipation in the external circuit cannot be neglected (C 0 *0). 


Table 1 

Numerical values of the model system. 


System parameter 

Symbol 

Value 

Resonator diameter 

Dr 

0.102 m 

Length ac resonator 

Lac 

2 m 

Regenerator diameter 

D r 

0.0889 m 

Length of regenerator 

L r 

0.073 m 

Specific impedance 

Zr 

3.6 x 10 9 m -2 

Length of pulse tube 

U 

0.24 m 

Diameter of pulse tube 

D t 

0.078 m 

Average length of space d 

L d0 

0.209 m 

Diameter of space d 

D d 

0.085 m 

Length of inertance tube 

C 

0.256 m 

Diameter of inertance tube 

Di 

0.078 m 

Average volume of space c 

V C 0 

0.00283 m 3 

Ambient temperature 

T a 

300 K 

Average pressure 

Po 

3 MPa 

Specific-heat ratio 

y 

1.67 

Viscosity at T a 

Pa 

20 ps Pa 

Density 

Po 

4.81 kg/m 3 

Conductance orifice 

Co 

O.lCr 
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4. Some characteristic dependences 

4.1. Input parameters 

In this section some examples of calculated results will be given which show what can be done with the results obtained 
so far. It is not the purpose of this section to give an exhaustive treatment of all implications of the formulae given above. 
Unfortunately Ref. [1] does not give enough detailed information to make a quantitative comparison between the results, 
derived in this paper, and the actual performance of the system. The numerical values of our model system are given in 
Table 1. They are the same as the system described in Ref. [1]. 

4.2. Some plots 

Next some plots are given for t c and the scaled oscillation frequency v/v ref with v ref = co eR /2n. They are plotted as 
functions of C 0 /C r , Dj, and L ac while the other system parameters have fixed values given in Table \. In Fig. 4 the 
conductance C 0 of the orifice is varied. We see that an onset temperature above ambient is necessary (t c > 1) even if C 0 = 0. 



0.0 0.1 0.2 
C/C 

o r 


Fig. 4. Onset reduced temperature z c and scaled oscillation frequency v/v ref as functions of the relative orifice conductance. 



0.00 0.05 0.10 

D i [m] 

Fig. 5. Onset reduced temperature t c and scaled oscillation frequency v/v re f as functions of the inertance diameter D,-. 



L [m] 

ac L J 

Fig. 6. Onset reduced temperature t c and scaled oscillation frequency v/v ref as functions of the resonator length L ac . 
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For stable oscillations a higher temperature is needed if the orifice conductance is increased. The oscillation frequency does 
not depend strongly on the orifice setting. 

Fig. 5 is a plot with the inertance diameter D, as variable. This plot shows a clear minimum of t c for D t = 4.5 cm. Fig. 6 is 
a plot of the onset temperature and the reduced oscillation frequency as functions of L ac . In this case there is a clear 
minimum of t c at L ac - 0.2 m. 


5. Transient behavior 


So far T t has been treated as an input parameter which is constant in time. We have seen that T t <T c corresponds with 
a<0, so oscillations tend to die out. On the other hand T t >T c corresponds with a>0. In this case the oscillations would 
grow without limitation if T t is fixed. Hence also the heat flow Qt, needed to fix T t , also would grow without limitation. It is 
more realistic to fix the heat flow Qt. This case will be considered in this section. 

In the presence of an oscillating gas flow in the engine the temperature T t oscillates with the same frequency. However, 
at frequencies in the acoustic region, the amplitude of the oscillations in the temperature of the heat exchanger is very 
small due to its high heat capacity C h . If only slow temperature variations are considered, the energy balance of the heat 
exchanger is given by 


d Tt • • * 

C H = Qt ~ Qc ~ H t . (51) 

Here Qt the applied heating power and Q c the heat flow conducted to the surroundings. If the only path of the heat flow is 
through the regenerator then 

Qc = K ay~(Tt — Ta) (52) 

with K a the effective thermal conductivity of the regenerator. The average enthalpy flow in the pulse tube is taken over one 
oscillation period so we write 




(53) 


where t p = 1 /v is the time between two zero crossings of 3p t with a positive slope. Defined in this way, the average 
enthalpy flow is a short-time average which can vary slowly with time. The pressure amplitude p 1 is defined as the 
maximum value of Sp t during the time span t p . 


With 


and 


Eqs. (6) and (12) result in 


and 


dV d 

d t Jd 

(54) 

Q- 

^ 5? 

ii 

• 

(55) 

d/d 

dt = a ' Pr 

(56) 

& = a « Sp <- 

(57) 


Eqs. (32), (34), (56), and (57) form a complete set of four first-order differential equations. The general solution of such a set 
is described in Ref. [6]. Here the set is solved numerically by second-order Runge-Kutta integration [7]. 

Fig. 7 gives the calculated time dependence of T t and the amplitude of the pressure oscillations in the pulse tube p 1 . 1 
The value of the thermal conductance of the regenerator K a A r /L r was set at 0.085 W/K. In order to speed up the effects the 
heat capacity C H was given the rather low value of 0.21 J/K. For this system the onset temperature T c = t c T a = 802 K. This is 
represented by a horizontal line in Fig. 7. With a thermal conductance of 0.085 W/I< a critical heating power of 43 W is 

needed to maintain a temperature of 802 K. The applied heating power was Qt = 500 W, so far above the critical value. 

* 

With this heating power the equilibrium temperature without oscillations (so with H t = 0) would be as high as 6182 K, so 
far above the onset temperature of 802 K. 

Now let us have a look at Fig. 7. At t = 0 the initial temperature T t was set at 6001< and Sp t was given a small kick of 
50hPa. As long as T t < 802 K the temperature of the hot end is too low to sustain steady oscillations (a<0). The pressure 


1 In this discussion we will only mention p^ as the amplitude of the oscillations, but it should be understood that the amplitudes of all other dynamic 
parameters, such as the pressure drop over the regenerator and the various volume flows, vary in the same way. 
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f[s] 

Fig. 7. Time dependences of T t (full line) and p^ (dotted line) for an artificially low C H . The horizontal line represents the critical temperature. 



Fig. 8. Time dependences of T t (full line) and p 1 (dotted line) for a realistic C H . The horizontal line represents the critical temperature. The figure shows a 
repetition of bursts of high-intensity sound interrupted by relatively long quiet periods. 


oscillations, resulting from the kick at the start, tend to die out. This is visible between t = 0 and 0.1 s. However, T t increases 
rapidly and after about 0.1 second T t reaches T c = 802 K. From then on T t >T c and the amplitude of the oscillations grows 
(a>0). This is associated with a rapid increase of the sound intensity illustrated in Fig. 7 by an increase of p^. The frequency 
v of the sound is about 103 Hz. Subsequently the hot end is cooled strongly by the oscillating gas flow (expressed in Eq. (51) 
by a high value of the average enthalpy flow) and, after going through a maximum, T t decreases. But the amplitude of the 
sound keeps on growing as long as T t >T c . At the moment that the hot end has cooled to the onset temperature T c the 
intensity of the oscillations is at its maximum! So T t continues to drop to values below T c , initially with a slow decay but 
later with a much stronger decay until the cooling of the hot heat exchanger by the oscillating gas becomes smaller than the 
applied heating power. From that moment the temperature T t starts to move up again. 

This is repeated with decreasing excursions until a state of steady oscillations is reached at T t = 802 K. This value, 
obtained from the numerical integration, corresponds with t c = 2.67 which is in excellent agreement with the value 
obtained from the analytical treatment (Eq. (47)). Also the numerically-obtained frequency v of 102.9 Hz is equal to the 

* 

value obtained by the analytical solution (45). In the steady state the enthalpy flow H t =459W and the heat flow by 
conduction is Q c = 41 W. 

The situation for T t « T c can be understood as follows: if T t <T c the oscillations are damped. Due to the applied heating 
power the temperature rises. At the moment T t passes T c the oscillations grow according to exp(at), but a grows linearly in 
time. So the oscillation amplitudes grow more than exponential. Setting t = 0 at the moment that T t = T C the oscillations 
grow with exp(aT 2 dT t /dt) with a' a positive constant. 

The repetition rate of bursts of acoustic power (with a frequency of about 2 Hz in the example of Fig. 7) is determined by 
the interplay between the heating power, the rate of the grow and decay of the amplitude of the high-frequency oscillations 
(sound intensity), and the heat capacity C h of the region around the hot heat exchanger. At values of the heating power just 
above the critical value the repetition rate is small and increases with Qt- 

Fig. 8 represents the calculated variations of T t and p^ as functions of time for the more realistic value of C H = 21J/K. 
The applied heating power was 2 kW and the starting temperature 750 K, so not too far below the critical temperature. The 
figure shows a repetition of bursts of very high intensity sound interrupted by relatively long quiet periods. Also in this case 
we see transient effects with a certain time period which is, in this case, about 3 s. The oscillation amplitudes of the 
pressure are strong functions of time and, in this example, can vary over an order of magnitude. 
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It is interesting to note that, under the conditions in this paper, the value of T c is determined only by the geometry of the 

machine and not by the heating power Qt . Consequently, in the steady state, also the heat flow Q c is independent of Qt- 

. * 

Only the oscillation amplitude p 1 is affected by Qt and, consequently, the average enthalpy flow H t . During the transient 
period, at the peak values of the p\ -oscillations, the energy flows are very high (in the example of Fig. 8 more than 12 kW), 
so large radial temperature gradients will be present heat exchanger. This will affect the dynamics of the behavior. Also the 
assumptions of linearity of the relations will not be valid at high amplitudes. It is an interesting situation that the onset of 
the oscillations is determined by the values of, e.g. C R and C 0 for small amplitudes, but that the effective values of these 
system parameters may differ for oscillations with high amplitudes. 

The transient behavior, reported in this section, has been observed experimentally [9]. Penelet et al. have reported 
effects that strongly resemble the pulsating oscillations described here [10]. Penelet et al. conclude that their model does 
not reproduce the experimental results very well. They argue that the most critical assumption in the model is the one¬ 
dimensional approach in the description of the heat transfer. It would be interesting to reanalyze their system in terms of 
the model presented in this paper. Yu et al. used CFD software to calculate the transient behavior in their machine for two 
conditions: one with fixed heat input (case A) and one with fixed hot-end temperature (case B) [11]. Case A corresponds 
with the numerical approach in this paper (Section 5). They find an overshoot of the pressure amplitude before it tends to a 
state of steady oscillations. This is resembles the results found here. Case B of Yu et al. corresponds with the analytical 
treatment of this paper. In the treatment, given in Section 4, the amplitude would grow without limitation once the critical 
temperature is passed. This is in contrast with the findings of Yu et al. The reason is most probably that the dissipation 
strongly increases at high amplitudes, e.g. due to vortex production. This is favored by the geometry chosen by Yu et al. 
where the gas is forced to make a 180° angle. This kind of effects is not taken into account in this paper (although is can be 
incorporated in C 0 ). Furthermore heat exchange may become problematic if the gas oscillates with high amplitude so the 
condition of constant temperature is not satisfied. The combination of these effects would also explain that there is only a 
single small overshoot in Case A. 

6. Discussion 

In this paper the basic features of thermoacoustic engines are treated with the travel ling-wave thermoacoustic heat 
engine as a working model. If the hot-end temperature is assumed to be constant the dynamics of the system is described 
by a fourth-order differential equation which determines damping, growth, or stability of oscillations. For the case that the 
hot-end temperature is allowed to vary in time the dynamic relations are integrated numerically. Interesting transient 
effects are seen in which the amplitude of the oscillations is pulsating. Eventually a steady state is reached in which the 
amplitudes of the oscillations are stable. The results of the numerical integration for the steady state agree very well with 
the analytical relations. 
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Appendix A. Derivations 

Eqs. (35) and (36) with the following identifications: 

a = (t t w e + w c )C r , b = (w c + w e )a u c = -w e C 0 , f = -a R w e , 

k = w e C 0 , l = a R w e , m = -i t w e C r , n -apN e 

become 


and 



Multiplying Eq. (61) with the operator 


d 2 dp t A5p t d Pr 

~dT Pt ~ ~df Pn 


n _ d d h 

a- df2 +a dt + fa 


(58) 

(59) 


(60) 


(61) 


( 62 ) 
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gives 


or 


With Eq. (60) 




+ k 


d Sft 
dt 


+ 





+ np r 










Expanding this relation and collecting terms gives 



with 


o _ d4s Pt , a d3s Pt , „ d 2 <5 Pt d dp t , 

U — —:—^-r U3 ——■5-r Q 2 — r 6q —~tt r ugOpf 


dt 4 


dt 3 


dt 2 


dt 


— k -t - a, 


(63) 

(64) 

(65) 

( 66 ) 

(67) 


a 2 = l + ak - cm + b, 


( 68 ) 


a 1 = al-fm - cn + bk , 


(69) 


a 0 = bl-fn. (70) 

If we would eliminate 3p t instead of p R we would have obtained the same differential equation but now in p R . This reflects 
the fact that the stability condition for all dynamic variables is the same. 

We write Eq. (66) as 


n _ d 4 <5p t a 3 d 3 <5p t a 2 d 2 Sp t eg dSp r 

a 0 dt 4 % dt 3 Oq dt 2 Oq dt t 

This suggests a characteristic frequency 


We introduce the dimensionless time x according to 


60 0 = 0 . 


1/4 

0 


then Eq. (71) reduces to 


with 


X = 60 0 1 


0 = 


d 4 <5p t d 3 <5p f d 2 Sp t 
dx 4 3 dx 3 2 dx 2 


+ C 1 



(71) 


(72) 

(73) 


(74) 


c 3 = 


CL 3 


c 2 = 


a 2 

60 ? 


Cl = 


Q .| 

CD} 


^0 ^0 

Note that all coefficients in Eq. (74) are real and positive. The solution of Eq. (74) is of the form [7] 


(75) 


4 

dp t = J2C k exp(z k x), (76) 

k= 1 

where the z k are the roots of the characteristic equation 

0 = z 4 + c 3 z 3 + c 2 z 2 + c-(Z+ 1. (77) 

Eq. (76) is the general solution of the time dependence of Sp t for constant r t . It contains not only eventual oscillatory 
solutions, but also all transient effects such as the transition to a steady state after an arbitrary starting condition. Eventual 
oscillatory solutions are not a prerequisite of our treatment, but rather the result. 

Eq. (77) can be solved analytically [8]. If one root is complex its complex conjugate is also a root since the coefficients , 
c 2 , and c 3 are real. In that case the linearly independent solutions are 

e ax cos 6 o x x and e ax sin co x x. 

The index x in co x is used to indicate that co x is the angular frequency expressed in the dimensionless time x. 


( 78 ) 
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We substitute 


z = tx + ico x (79) 

in Eq. (77). We are interested in the onset of oscillations, so we can limit the discussion to low values of a. To first order in a 
Eq. (77) gives 


0 — -4iaco 3 - 3c 3 acOx + 2ic 2 aco x + CjOC + co# - ic 3 co 3 - c 2 co^ + iCjCO* + 1. 

The real part of Eq. (80) results in 


(80) 


and the imaginary part in 


0 = cd x - (3c 3 a + c 2 )cl>x + Cj a + 1 


Substituting Eq. (82) in (81) gives, with Eq. (75), to lowest order in a, 


a = 


1 C 1 C 3_!_ M 03 

2 c\ + (2c-[ - c 2 c 3 ) 2 of 2 \0 3 0 Qi 



As the prefactor is positive the sign of a is determined by the sign of the function 


s = 


cq 

03 


+ a 0 


0 3 

cq 


a 2 . 


If s<0 the oscillations are damped. If s>0 the oscillation amplitude is growing. 
For steady oscillations s = 0. In this case Eq. (82) gives 


The real-time angular frequency is given by 



c 3 


With a = 0 Eqs. (81) and (82) give 


Substitution of Eq. (87) in Eq. (77) reads 



2 2 

= COqCOx 


cq 

«3* 


Cl c 3 
c 2 = — + — 
c 3 cq 


Eq. (88) can be written as 


0 = z 4 + c 3 z 3 + 





+ Ci z + 1. 


°=( z 2 + S)( z2+C3Z+ ^)- 

The first factor gives the two roots 



which we know from Eq. (85). The second set of roots z 3 4 is given by 




(81) 

(82) 

(83) 

(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 

(90) 

(91) 


As c 3 >0 the second set of solutions z 3 4 corresponds with damped oscillations. So there is only one steady oscillation 
possible. The property differs from real thermoacoustic systems where higher harmonics are possible. This can be 
mimicked in our model by choosing different values for the parameter a R . 

Table 2 gives the calculated coefficients in the characteristic equations of the model system. The fact that the 
dimensionless parameters q, c 2 , and c 3 are of order 1 is an indication that they are good dimensionless parameters to 
characterize the basic properties of the system. 
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Table 2 





Coefficients in the characteristic equations of the model system. 




Symbol 

Expression 

Value 

Symbol 

Value 

a 0 

71 

3.59 x 10 12 s -4 

c 0 

1 

a i 

70 

1.41 x 10 9 s“ 3 

Cl 

0.54 


69 

8.99 x 10 6 s“ 2 

C2 

4.75 

a 3 

68 

3.37 x 10 3 s- 1 

c 3 

2.45 

Table 3 





Characteristic frequencies for the model system defined in Table 1. 




Symbol 


Expression 


Value (Hz) 

v Ccr 


w c C R /2n 


332 

v Cer 


w e C R /2n 


73.8 

v Ceo 


w e C 0 /2n 


7.38 

v eR 


^a^W7/2n 


115 

Vci 


^0^/271 


417 

v ei 


^a(W7/2n 


196 

V 


Eq.(47) 


102.92 

Appendix B. 

Characteristic frequencies 




The combination of the compliance volume V c and the regenerator conductance C r , 

the volume V e and C r , 

and V e and 

the orifice conductance C 0 are all examples of case c of Fig. 2. The rate of change of the pressure in the vessel is determined 

by three characteristic frequencies 






&>Ccr — , 


(92) 



02cer = 


(93) 



(!) Ceo = w eQ)- 


(94) 


The product w e a R 


col R = a R w e = a R ^p- (95) 

v e 

is the square of the angular resonance frequency of the mass of the piston and the gas-spring action due to the volume V e . 
Similarly 


(D 2 d = OjWc, 


(96) 


CD 2 ei = OjWe (97) 

are the resonance frequencies of the mass-spring systems formed by the inertance and the compliance and the inertance 
and V e , respectively. Table 3 gives the values of the various frequencies for the model system. 

With these frequencies Eqs. (35) and (36) become 


d 2 p r 

dt 2 


+ i^t^Cer "T ^Ccr) 


dPr 

dt 


+(«;?,•+cu2.)p r 


d<5p t 

-®Ceo V - 



(98) 


and 


d 2 Sp t 

dt 2 


+ °^Ceo 


d Spt 
dt 


+ CD 2 R Sp t 


^t^Cer 


dp, 

dt 



(99) 


The treatment in this paper could be completely formulated in terms of the characteristic frequencies, but the relations are 
easier to recognize if they are expressed in terms of the conductances C r and C 0 , the masses M, and Mr, etc. which makes 
the paper easier to read. 
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